clear all;
clc;

randaxo = randperm(noaxons); % random permutation 

% ---------------------- Import Cord Geometry -----------------------------

patnum = 3;
cordgeom = load('patient_cordgeom.txt');

ae = cordgeom(patnum,1); % transverse left-right axis of white matter ellipse
be = cordgeom(patnum,2); % transverse down-up axis of white matter ellipse
c_white = [0;cordgeom(patnum,3);0]; % center of white matter ellipse

dr_th1 = cordgeom(patnum,4);
dr_th2 = cordgeom(patnum,5);
dc_th1 = dr_th2;
dc_th2 = cordgeom(patnum,6);

yf = cordgeom(patnum,7);

% XY coordinates
droot_xy = load(['p',num2str(patnum),'_dr_xyloc.txt']);
dcolm_xy = load(['p',num2str(patnum),'_dc_xyloc.txt']);


delr_dura = 0.2;
patientgeom = load('femspine_patientgeom.txt');
betageom = patientgeom(patnum,:);
a_cord = betageom(1)/2;
b_cord = betageom(2)/2;
a_csf  = betageom(3)/2-delr_dura;
b_csf  = betageom(4)/2-delr_dura;
a_dura = a_csf + delr_dura;
b_dura = b_csf + delr_dura;
ycenter = 13.61;
csf_center = [0,ycenter];
ytopepi = 25;

y_cord = (ycenter-b_csf) + betageom(5)*2*b_csf;
cord_center = [0,y_cord];

% ---------------------------- Boundaries ---------------------------------

% spinal cord boundary
s  = 0:pi/100:2*pi;
r_ellipse = ae*be./sqrt((be*cos(s)).^2+(ae*sin(s)).^2);
x_ellipse = r_ellipse.*cos(s);
y_ellipse = r_ellipse.*sin(s);

% grey matter
xg = [0.00 ,0.62 ,1.24 ,1.87 ,2.09 ,2.24 ,1.45 ,2.03 ,2.22 ,2.03 ,1.21 ,0.62 ,0.00];
yg = [-0.44,-0.12,0.96,2.04,2.06,1.88,0.43,-1.12,-1.62,-2.12,-2.18,-1.61,-0.88];
yg = yg;
xgrey = [xg,-fliplr(xg)];
ygrey = [yg,fliplr(yg)];

% dura boundary
r_dura= a_dura*b_dura./sqrt((b_dura*cos(s)).^2+(a_dura*sin(s)).^2);
x_dura = r_dura.*cos(s);
y_dura = r_dura.*sin(s);

% -------------------- Electrode Locations --------------------------------

% lateral (clockwise) deviations of 0, 10, and 20 degrees
% - converted to radians
th1 = (90-0)*pi/180;
th2 = (90-10)*pi/180;
th3 = (90-20)*pi/180;

% 3 epidural locations
aepi = a_dura + 1;
bepi = b_dura + 1;

repi1 = aepi*bepi./sqrt((bepi*cos(th1))^2+(aepi*sin(th1))^2);
xepi1 = repi1*cos(th1);
yepi1 = repi1*sin(th1);

repi2 = aepi*bepi./sqrt((bepi*cos(th2))^2+(aepi*sin(th2))^2);
xepi2 = repi2*cos(th2);
yepi2 = repi2*sin(th2);

repi3 = aepi*bepi./sqrt((bepi*cos(th3))^2+(aepi*sin(th3))^2);
xepi3 = repi3*cos(th3);
yepi3 = repi3*sin(th3);

xepi = [xepi1,xepi2,xepi3];
yepi = [yepi1,yepi2,yepi3];

% 6 subdural locations

% 1 - 1 mm above cord
asub1 = a_cord + 1;
bsub1 = b_cord + 1;

rsub11 = asub1*bsub1./sqrt((bsub1*cos(th1))^2+(asub1*sin(th1))^2);
xsub11 = rsub11*cos(th1);
ysub11 = rsub11*sin(th1);

rsub12 = asub1*bsub1./sqrt((bsub1*cos(th2))^2+(asub1*sin(th2))^2);
xsub12 = rsub12*cos(th2);
ysub12 = rsub12*sin(th2);

rsub13 = asub1*bsub1./sqrt((bsub1*cos(th3))^2+(asub1*sin(th3))^2);
xsub13 = rsub13*cos(th3);
ysub13 = rsub13*sin(th3);

% 2 - 1 mm below dura
asub2 = a_dura - 1;
bsub2 = b_dura - 1;

rsub21 = asub2*bsub2./sqrt((bsub2*cos(th1))^2+(asub2*sin(th1))^2);
xsub21 = rsub21*cos(th1);
ysub21 = rsub21*sin(th1);

rsub22 = asub2*bsub2./sqrt((bsub2*cos(th2))^2+(asub2*sin(th2))^2);
xsub22 = rsub22*cos(th2);
ysub22 = rsub22*sin(th2);

rsub23 = asub2*bsub2./sqrt((bsub2*cos(th3))^2+(asub2*sin(th3))^2);
xsub23 = rsub23*cos(th3);
ysub23 = rsub23*sin(th3);

xsub = [xsub11,xsub12,xsub13,xsub21,xsub22,xsub23];
ysub = [ysub11,ysub12,ysub13,ysub21,ysub22,ysub23];

% -------------------- Transverse DC/DR Seeding ---------------------------

noaxons = 200;

% DC parameters
nodclv    = 10;
axonperlv = 10;
nodiv     = nodclv;
ptsperdiv = axonperlv;
nopts = nodiv*ptsperdiv;

% seed_dclayer(nodclv,axonperlv,dc_th1,dc_th2,0,yf*be,ae,be);
% seed_dclayer(nodiv,ptsperdiv,tho,thf,yo,yf,ae,be)

% convert to radians
tho = dc_th1*pi/180;
thf = dc_th2*pi/180;
ybot = 0;
ytop = yf*be;

seed_coord = zeros(nopts,2); % 2 columns: x,y

% splitting arc azimuthally
th_topbnd = tho:(thf-tho)/nodiv:thf;
rtopbnd = ae*be./sqrt((be*cos(th_topbnd)).^2+(ae*sin(th_topbnd)).^2);
xtopbnd = rtopbnd.*cos(th_topbnd);
ytopbnd = rtopbnd.*sin(th_topbnd);
% splitting down-up radius of ellipse
ybnd = ybot:(ytop-ybot)/nodiv:ytop;

% -------------------- Plotting Transverse Seeding ------------------------

figure;
hold on;
plot(x_dura,y_dura,'k--','LineWidth',4);
fill(x_ellipse,y_ellipse,[1,1,1],'LineWidth',3);
fill(xgrey,ygrey,[0.5,0.5,0.5],'LineWidth',3);
plot(xepi,yepi,'ks','MarkerSize',12,'LineWidth',3);
plot(xsub,ysub,'k^','MarkerSize',12,'LineWidth',3);
% (100,100) are dummy points, not on plot, for the sake of legend
plot(100,100,'k.','MarkerSize',30);
plot(100,100,'kd','MarkerSize',12,'MarkerFaceColor','k');
hold off;
L1 = 'Dura Mater';
L2 = 'White Matter';
L3 = 'Grey Matter';
L4 = 'Epidural Locations';
L5 = 'Intradural Locations';
L6 = 'DC Fibers';
L7 = 'DR Fibers';
legend(L1,L2,L3,L4,L5,L6,L7,'Location','NE');
xlabel('x (mm)','FontSize',30);
ylabel('y (mm)','FontSize',30);
set(gca,'FontSize',30);
axis([-10,10,-10,10]);
axis square;

figure;
hold on;
fill(x_ellipse,y_ellipse,[1,1,1],'LineWidth',3);
fill(xgrey,ygrey,[0.5,0.5,0.5],'LineWidth',3);
for ii = 1:(nodiv+1)
    x1 = 0;
    y1 = ybnd(ii);
    x2 = xtopbnd(ii);
    y2 = ytopbnd(ii);
    plot([-x2,x1,x2],[y2,y1,y2],'k-','LineWidth',3);
end
sub = 1:7:noaxons;
plot(dcolm_xy(:,1),dcolm_xy(:,2),'k.','MarkerSize',20);
plot(droot_xy(sub,1),droot_xy(sub,2),'kd','MarkerSize',10,'MarkerFaceColor','k');

xlabel('x (mm)','FontSize',30);
ylabel('y (mm)','FontSize',30);
axis([-2,2,-1,3]);
axis square;
set(gca,'FontSize',30);

% ------------------------ Fiber Data -------------------------------------

droot_data = load('p1_drootxyz_200axons_167nodes_9um.txt');
dcolm_data = load('p1_dcolmxyz_200axons_312nodes_9um.txt');
droot_data(2,:) = droot_data(2,:)-y_cord;
dcolm_data(2,:) = dcolm_data(2,:)-y_cord;

elepernode = 8;
eleperdc = (312-1)*elepernode+1;
eleperdr = (167-1)*elepernode+1;

% split dorsal column x, y, and z
dc_x=zeros(eleperdc,noaxons);
dc_x(:)=dcolm_data(1,:);
dc_x=dc_x';
dc_y=zeros(eleperdc,noaxons);
dc_y(:)=dcolm_data(2,:);
dc_y=dc_y';
dc_z=zeros(eleperdc,noaxons);
dc_z(:)=dcolm_data(3,:);
dc_z=dc_z';

% split dorsal root x, y, and z
dr_x=zeros(eleperdr,noaxons);
dr_x(:)=droot_data(1,:);
dr_x=dr_x';
dr_y=zeros(eleperdr,noaxons);
dr_y(:)=droot_data(2,:);
dr_y=dr_y';
dr_z=zeros(eleperdr,noaxons);
dr_z(:)=droot_data(3,:);
dr_z=dr_z';

% ---------------------- Truncating Coordinates ---------------------------

ss = 25;

dc_coord=cell(noaxons,1);
dr_coord=cell(noaxons,1);
for ii=1:noaxons
    
    ck_dc=( abs(dc_z(ii,:)) <= ss );
    dc_coord{ii}=[dc_x(ii,ck_dc);dc_y(ii,ck_dc);dc_z(ii,ck_dc)];
    
    ck_dr=( abs(dr_z(ii,:)) <= ss );
    dr_coord{ii}=[dr_x(ii,ck_dr);dr_y(ii,ck_dr);dr_z(ii,ck_dr)];
    
end


% ---------------------- 3D Plot of DC/DR Fibers --------------------------

[xc,yc,zc] = cylinder([1,1],20);
xc = (ae+0.2)*xc;
yc = (be+0.2)*yc;
zc = ss*(2*zc-1);

xtmp=xc;
ytmp=yc;
ztmp=zc;
xc=xtmp;
yc=ztmp;
zc=ytmp;

figure;
hold on;
surf(xc,yc,zc,'FaceColor',[0.5,0.5,0.5],'EdgeColor','none');
colormap hsv;
alpha(0.3);
grid off;


subsamp_dc = randaxo([1,20:20:noaxons]);
% subsamp_dr = randaxo([1,5:5:noaxons]);
% subsamp_dc = [1,5:10:noaxons];
subsamp_dr = [1,5:5:noaxons];
for ii = subsamp_dc
    plot3(dc_coord{ii}(1,:),dc_coord{ii}(3,:),dc_coord{ii}(2,:),...
        'k--','LineWidth',2);  
end
% for ii = subsamp_dr
%     plot3(dr_coord{ii}(1,:),dr_coord{ii}(3,:),dr_coord{ii}(2,:),...
%         'r-','LineWidth',2.5);    
% end
% fill3(xgrey,ss+zeros(size(xgrey)),ygrey,[0.5,0.5,0.5],'LineWidth',1.5);
hold off;
axis([-(ss+2.5),ss+2.5,-(ss+2.5),ss+2.5,-(ss+2.5),ss+2.5]);
axis on;
axis square;
axis off;
% view([-90,0]);
view([0,90]);

xlabel('x (mm)','FontSize',30);
ylabel('y (mm)','FontSize',30);
zlabel('z (mm)','FontSize',30);
legend('Spinal Cord','DC Fiber','DR Fiber','Location','NE');
set(gca,'FontSize',30);

